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Abstract 


The behavior of pressure transient during production tests in radially composite multilayered reservoirs 
with non communicating layers has been widely studied, as well as multilayered reservoirs with homogeneous 
layers that communicate in the formation. In that sense, the present work attempts to develop an analytical 
model to calculate the pressure response for multilayered radially composite reservoirs considering the 
presence of formation crossflow along all layers. 

An infinite reservoir will be considered and the properties for each layer, such as permeability, oil viscosity, 
and porosity may be different. The analytical model proposed in this work, was developed based on 
equations that model pressure behavior during conventional production tests in multilayered reservoirs. 
Moreover, the mathematical formulation for this model is presented in the Laplace domain and the Stehfest 
Algorithm is used to convert the pressure variation solution into the real field. The results presented in this 
work showed a close agreement when compared to the responses provided by commercial flow simulator 
based on finite differences. 


1. Introduction 


There are different models of oil reservoir to consider during a production test. In this work, the system 
considered is a multilayered reservoir with different regions of permeability in each layer and, in addition, 
formation crossflow is considered. Hence, there is a nonzero vertical permeability between two adjacent 
layers. Formation crossflow is different than crossflow in the wellbore, which is also considered. Formation 
crossflow is defined in this work as in the semipermeable-wall system of Sun and Gao (2017). 

One important work regarding single-phase flow in stratified reservoir system, which inspired many 
others, was the one by Lefkovits et al. (1961). In this work, an analytical model was proposed, where the 
properties, such as thickness, permeability, and porosity could be different in each layer. Formation damage 
was also considered. However, a commingled system was considered, i.e, formation crossflow is not included. 
On the other hand, the works presented by Ehlig-Economides and Joseph (1987), Bourdet et al. (1985), 
Prijambodo et al. (1985) and Sun and Gao (2017) included formation crossflow. 

The work proposed by Bourdet et al. (1985), was entirely done in the real field. This work proposes 
an analytical model considering formation crossflow, for a two layer reservoir. It also considers wellbore 
storage and skin effects. In addition, it includes a double porosity and permeability model description. 

The work proposed by Ehlig-Economides and Joseph (1987) was strongly based on the previous Lefkovits 
et al. (1961) and Sun and Gao (2017). In this work, an analytical model is presented considering a reservoir 
with an arbitrary number of layers under single phase flow, however this work included formation crossflow. 
Properties may also be different from layer to layer. The presented model considers both formation damage 
and wellbore storage, and a derivation was made for both short and long times. This work showed that after 
a period of time, the reservoir could be described as an equivalent single layer system. This very complete 
work was the main study source for this present work. The formulation regarding the interflow between 
layers was derived from it. In addition, the work presented by Ehlig-Economides and Joseph (1987), unlike 
Bourdet et al. (1985), was developed in Laplace field and used numerical inversion of Laplace transforms, 
more specifically, the Stehfest Algorithm (Stehfest, 1970) in order to find a pressure response in the real 
field. The first work to use it, was the one proposed by Tariq et al. (1978). 

The study presented by Prijambodo et al. (1985), considered a two layer cylindrical reservoir. It also 
considered damaged regions for the results. It presents an analysis for producer wells. It showed that the 
flowing pressure response at the well at early times can be divided in three flow periods; The first one shows 
that the reservoir behaves as if there were no crossflow, the second one, a transitional period, shows that the 
pressure response depends on the contrast in horizontal permeabilities and on the degree of communication 
between layers. During the third period, the reservoir behaves like an equivalent to a single-layer system, 
like the results presented in Ehlig-Economides and Joseph (1987) and Sun and Gao (2017). 

Moreover, a radially composite reservoir is considered in this work, i.e, the presence of different regions 
of permeability within each layer. The work of Closmann et al. (1967) considered a composite radial single 
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layered reservoir’s model. Likewise, the work proposed by LF Bittencourt Neto et al. (2020) considered 
a composite radial reservoir, however a two-layered model was considered. The formulation done in 
(LF Bittencourt Neto et al., 2020) was, like in this work, developed in the Laplace domain. Both of these 
works (LF Bittencourt Neto et al., 2020) and (Closmann et al., 1967) did not considered crossflow. 


There has been many studies regarding multiple regions and transient time and also, many studies 
regarding the behavior of a multilayered reservoir with formation crossflow alone as referenced before. 
However, this work is innovating as it combines both multiple regions within each layer and formation 
crossflow along the layers’ adjacencies. In addition, the model proposed can be used to calculate the 
equivalent permeability of the reservoir and this model was validated by comparison with a commercial 
flow simulator. 


The next section of this paper, section 2, describes the model’s hypothesis. Section 3 presents the 
mathematical analysis of the problem. Section 4 is reserved for results and discussions, and finally, section 
5, presents the conclusions. 


2. Model Description 


Consider an infinite, laterally isotropic, multilayered, radially composed reservoir. In each one of its n 
layers denoted by 7, there are m; distinct regions of permeability denoted by 7, and there may be formation 
crossflow in between all layers’ adjacencies. All the formulation was done in consistent units. 


Besides the permeability (k;;), each layer is considered homogeneous to all other properties. It has a single 
phase flow of a fluid with viscosity (1) that is constant, and also constant and small total compressibility 
(c,) in all layers. In addition, the porosity (¢;) may have a different value in each layer, as well as the 
thickness (h;). Constant flow rate at the well is considered. Formation damage and wellbore storage will be 
disregarded. Consider Apj; = pj — pji, the flow in each region 7 of each layer 7 = 1,2,...,n is governed by 
the following diffusion equation (Bourdet et al., 1985; Ehlig-Economides and Joseph, 1987): 


OAD ji 
Ot 





(kh) V7 Apyi = (Gh) jer + (Apyi — Apg—1i) X51 — (Apgar — Avy) X} [1] 


In order to define the coefficient of semipermeability X}, first consider that the radii of the regions of 
permeability are all equal for each layer j7. Figure 1 illustrates that problem for a 2 layers and 2 regions 
case: 
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dt 
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Fig. 1 — Two layered model with regions of equal radii 


Then, each coefficient of semipermeability Xi is defined as (Ehlig-Economides and Joseph, 1987): 





. 2 
XO BI(hy) hol] +38 + Xan e 
For j =1,...,n andi=1,...,m,. In this case, mj = mj41 for all j =1,...,n—1. 
Ah; and Key are the thickness and the permeability of the shale between layers j and j + 1, xi — ee 
where ki is the vertical permeability of layer 7 in region 7. ‘ 
Notice that Xj, = 0 and if there is no crossflow between layers j and j + 1, then X} = 0. 
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hy Layer 2 








Fig. 2 — Two layered model with regions of different radii 


Now, if the radii of adjacent regions in layers 7 and 7 + 1 are different, there will be coefficients of 
semipermeability that relates different regions. Figure 2 illustrates that problem for a case considering two 


layers and two regions. 
Notice that, in figure 1 the semipermeability coefficient X? relates region 2 from layer 1 and region 2 


from layer 2. Then, in that case, X? is given by: 


2 
Xe= 3 
1 2[Ahy/k2,] + x3 + x? 3) 





where y3 = ae and y? = Hh. 
22 ZY) 
However, in the example illustrated in figure 2, X? relates different regions of permeability. In this case, 
X? relates region 2 in layer 1 and region 1 in layer 2. Then, X? is given by: 





y: 
Xe= 4 
1 2[Ahy — /k2,.J+xh+ x? 4 


where x3 = a and x7 = yh. 
22 ZY 
Then, Xj} can be defined generally (for equal and different radii cases) as: 


. 2 
Xj = j rin [5] 
2[(Ah;)/kv5] ep eed 





where regions 7; and i;41 can represent equivalent or different regions for layers 7 and j + 1. 
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The example described in figure 2, considers three coefficients of permeability. Then, each interface can 
be extended to the other layer and the problem can be treated like one with three distinct regions in each 
layer, where the radii of the three regions are the same in both layers. 

Consider now, m; to be the number of regions in layer 7. From the examples discussed above, it is 
possible to see that even if the numbers of the regions are the same for two connected layers, i.e, m; = 
mj41, the radius of each one of these m,; regions can still be different. 

Let r(j, 1), r(j, 2), r(g,3),...,r(j,mj;) be the radii of the regions in layer j and for layer j + 1, r(j + 
Lae ly 2), Fj 18) s0cs PU AL a): 

Let s be the number of incidents where r(j,/) = r(j + 1,k), then, it is possible to show that the number 
of Xi semipermeability coefficients in the adjacency between layers 7 and j + 1 is m; +mj41 —1—s. 

For example, in figure 3, the first two layers represented have, respectively, three and four regions. Also, 
r(1,2) and r(2,3) are the only values that coincide (r(1,2) = r(2,3)). That is, m1 = 3, mz =4 and s=1. 
Hence, there are 3+ 4 — 1 —1=5 crossflow coefficients within these layers. 
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Layer n-1 
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| [ | | | | | Layern 
—\ —_. —=— 
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Fig. 3 — Reservoir with n connected layers and many regions of permeability 


After having each semipermeability coefficient well defined, consider the following equation: 


a= (Dim -1))-s [6] 
= 


Then, extending each interface to all reservoir’s layers and creating artificial regions, illustrated in figure 
4, this problem can be reduced to one with a regions of equal radii. 
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Fig. 4 — a extensions of the regions in each layer 


It is important to define well each semipermeability coefficient Xi for the a regions, i.e, observing 
correctly the permeability regions it is relating. 

After all layers have been divided into a@ regions, there might be regions with the same X 3 coefficients. 
For example, it is possible to see that in figure 4 X{ = X?, because both of the coefficients relate region 1 
in layer 1 with region 1 in layer 2. 

Notice that there are mj + mj41 — 1—s different values of X} and a — (mj +mj41 — 1 —) repeated 


values for each adjacency. 





3. Mathematical Formulation 


Considering the previous hypotheses and denoting m = a (number of regions), the model description is 
given by the diffusion equation [1] and its suitable boundary conditions: 


e The Initial Condition (I.C.) which occurs at time t = 0 and reflects that the reservoir is in balance, 
that is, the pressure is the same in all layers (apart from the hydrostatic effect). 


e The Internal Boundary Condition (I.B.C) which is related to the way in which the well is producing 
during the test. 


e The External Boundary Condition (E.B.C) which refers to the flow behavior at the extreme limit of 
the reservoir and in this work, an infinite-like reservoir is considered. 


Then, the following equation model the problem for all layers j = 2,...,n—1l and regions i = 1,...,(m—1): 
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OAp ‘4 i i 
PDE: (kh)jiV?Apyi = (bh) jicep— + (Apsi — Apg—1i) X51 — (Apgais — Adyi) X5 [7] 


When 7 = 1, the PDE is: 














OAD ;i 
Kyi Apyi = wy Ce a ~ = (Apgtiys — Aji) X5 [8] 
And for 7 = n, the PDE is: 
OAD ;; 
KyiV Apyi = Wy Crp os + (Ape = Apg ai) AG 4 [9] 
For all layers j = 1,...,n and all regions i = 1,...,(m— 1) there is the initial condition: 
IC: Ap;i(r,0) =0 [10] 
For all layers 7 = 1,...,n and region 7 = 1 there is an additional condition, the inner boundary condition 
(Lefkovits et al., 1961): 
. _ 2n(kh) jr OAD ji 
IBC: Di = ( Ap ) a [11] 
And for all layers 7 = 1,...,n and region 7 = m there is also an additional condition, the external 
boundary condition: 
EBC: im Apji(r,t) =0 [12] 
Consider 7 = 1,...,n and i = 2,...,m. There are coupling conditions relating the regions i — 1 and i 


in each layer 7, which must be defined in order to properly solve the problem. The coupling conditions 
between the regions (CCR) are given by the pressure and flow rate equality at the interface between them 
(Nie et al., 2011): 


Apj(i—1) (7 = 754-1), #) = Apgilr = 73-1), €) [13] 
Qj(i-1) = Qi 
Using Darcy’s law it is possible to rewrite the flow rate relation of the CCR so that : 
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[14] 





| Apji-1)(" = T34@-1)5 t) = Apgilr = 754-1); ) 


p OPI G-1) | — _tkh)ji_ (,, 0A; | 
Or T=Tj(i-1) a (kh) 546-1) or P= 5j(G—1) 


where j= 1,...,nandi=1l,...,m. 


The coupling condition between layers (CCL) is obtained considering that the pressure drop is equal in 
the layers along the well (except for the hydrostatic effect) and that the flow rate at the well is given by the 
sum of the flows of each layer, that is (LF Bittencourt Neto et al., 2020): 








Apg-11(r = Tw, t) = Apji(r = Tay b) [15] 
q= 1+ Gai +... + dni 
Once again, using the Darcy’s law, the flow rate CCL equation can be rewritten: 
Apgar — Popst) = Apyi(r = Tags't) 
qe x _ OAD; [16] 





2mrw Kyi Or |lr=ry 


In order to find the pressure solution, the Laplace transform will be applied in equations [7] through [16] 
and the solution will be first calculated in such space. Then, the Stehfest Algorithm (Stehfest, 1970) is 
used to find the solution in the real field. 

Consider Kj; = kjishj and w; = jh; to facilitate the analysis, Ap to be the Laplace transform of the 
pressure variation and u to be the Laplace’s variable. Now, for layers 7 = 1,...,n and regions 7 = 2,...,m-—1 
the following equation is given in the Laplace field: 








ODE: KyiV AD ji = wjycupAp;, + (AD;; — AD Gj—1)i) X}-1 — (Aya — AD,i)X} [17] 
For layers j = 1,...,n and region i = 1: 
_ Qnn, ru ,OAD,; 
IBC: qi = Zs ( Z) [18] 
Le OF ae, 
For layers j = 1,...,n and region i = m: 
EBC: lim Ap,;(7,t) = 0 [19] 


The coupling conditions between regions in the Laplace field are given below (LF Bittencourt Neto et al., 
2020; Nie et al., 2011): 
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AB, i_ay(t = TyG—1)>t) = ABT = ryi-1)st) 
H(i (" = P9194) = AP = P51 
CCR: ( ore) Kyi ( | 

rT — Or TST 5(4-1) — K5(i-1) r or r=T5(i-1) 


For j7=1,...,n andi=2,...,m. 


The coupling conditions between layers in the Laplace field are given below (LF Bittencourt Neto et al., 
2020): 


AD G—1(" = Tw, t) = Adi (7 = Tw, t) 
CCbey .— ts az (2922) 7 1] 
27U > Kj1 j=l r=Pw 
j=l 





For j = 2,...,7. 
The pressure general solution in the Laplace field is well known in terms of Bessel modified functions. 
They are given below for region i in layer j (Ehlig-Economides and Joseph, 1987): 


AD, = A‘ Ko(ojir) + Bi Io(oyir) [22] 
Applying the EBC in the equation above, for i = m: 
Jim [A Ko(ojer) + Bilo(ojir)] =0 > By =0 [23] 


Hence, the following equations for the general pressure solution are given: 


For regions 27 = 1,...,m—1: 


AD ji = A‘ Ko(ojir) + Bi Io(ojir) [24] 


For region 2 = m: 
AD; = Ai Ko(ojir) [25] 
Using the solutions above, it can be seen that V? AD; = oF, ADji- Hence, replacing this equivalence in 
the ODE [17] (Ehlig-Economides and Joseph, 1987): 
ji 5, ND ji — wiCeHudp,;, + (ADs; — ADGj—1)i) X}-1 = (Appia: = AD) X} =) [26] 
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And rearranging the equation above for each region i: 


(Kji05; — wyicywu — X}_4 — X))Apy + X}_ADG_1yi — XJ AD G41); = 0 [27] 


This is a homogeneous linear system for each region i where the nontrivial solution is wanted, that is, 
Ap,; # 0. That is true only if each matrix below is singular, and that implies that its determinant must 
vanish: 


aie fork=j—1;j >1, 
2 4 a 2 
Kjo° —wycpu—X; ,—-X:, fork= J, 
re i i — For i=1,...,m. (28) 
Aas fork =jtlsj<n, 
0, fork Aj—1,7,0or 7 +1. 


To find the values of oj; by vanishing the determinant of each matrix above is equivalent to finding the 
eigenvalues Kyio of each matrix below: 


—Xi4, fork =j—1;7>1, 
; wicput X?_,+X', fork =j, 
a a i jt : For i=1,...,m. [29] 
=KG) fork=jtlsj<n, 
0, fork Aj —1,7,0or 7 +1. 
This equivalence is true because: 
aig, = Dig — eig(biy,) I [30] 
where J is the identity matrix, and, consequently: 
det(a',) =0 = det( “ke — eig(bi,)I) =0 [31] 


Then, finding o;; is equivalent to extracting the square root of the eigenvalue of bi, then, dividing it by 
ae 


eig(bi,) 
Oxi = ne [32] 
Kai 
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Now, with the values of oj; calculated, the pressure coefficients can be determined. It can be seen in 
equations [24] and [25], that there are (2m — 1)n coefficients. The relations provided by the coupling 
conditions between layers and regions are, precisely, that many. 


The first n relations are provided by using the pressure solutions for region 1 [24] in the CCL pressure 
and flow rate equations: 


Pressure CCL: 
Aj_1Ko(oG-11Tw) + By_ylo(oG—11Tw) = Aj Ko(oj1Tw) + By Io(oj1T w) [33] 
For j = 2,...,n. 
Flow rate CCL: 


gy 
27Ury 





= oy1k11(—A}Ko(ourw) + Bilo(ourw)) +... + Onikni(—AlKo(onirw) + Blo(ontTw)) [34] 


The other (2m —1)n —n equations left are provided by using pressure solutions [24] and [25] in the CCR 
pressure and flow rate equations: 


Pressure CCR: 


Aj *Ko(o4@-1)"j(é-1)) a By *Io(o 40-1)" 5¢-1)) — Aj Ko( oi" j(i-1)) + Bylo(oiri¢-1)) [35] 


For 7 = 2,...,m—1. 


Aj, *Ko(o 4-1" i6-1) as By" To(o¢-1)36-1)) — Aj Ko(o5erj(i-1)) | 
For i= ™. 
Flow rate CCR: 
KjiO 54 ; i 
jeg (Ai Ka (opirj-1)) = Bil, (ojer@-1))) 
[37] 


Ak 0 5(4G-1) 7 (G— = Bray 05 (G-1) "7 (G— SS 
§ Ki(oya—1"3@-1)) — By Le j@—1"5@-1)) Pa Pa 


For 7 = 2,...,m—1. 


KyiO ji 


Ai Ki (056-17 96-1)) — BF (04-1) '3G-1)) = (Ap Ki (opirsi-1)) [38] 


Kj (i—1)9 j(é—-1) 
For i =m. 
Using the relations above, a linear system can be set up in order to calculate the pressure coefficients: 
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Aj qu/(2rrwu) 
Bi 


1 0 
“ 0 
By 0 

Ay" 0 
A 0 
B 0 
‘Be —mu- 0 
AS” 0 
An 0 
ay 0 
Ay 0 
a. 0 
Amn 0 


where the matrix M above is defined by the (2m — 1)n equations [35] to [38]. Finally, it is possible to 
find a solution in the Laplace field for the pressure variation at the wellbore: 


AD wf = A} Ko(ourw) + BlIo(ourw) [39] 


The Stehfest Algorithm (Stehfest, 1970) is used to invert this solution into the real field. 


4. Results and Discussions 


The results for the wellbore pressure profile for the proposed analytical model is presented in this 
section. The accuracy of the proposed solution was testified by comparing it with a commercial finite 
difference-based flow simulator (CMG, 2010). 

As an input for the numerical simulation, a radial grid was considered. The grid is more refined in the 
region closest to the wellbore, which is the region most affected by its presence. The oil model used in the 
simulator was blackoil. For all cases, formation crossflow is considered between the adjacent layers. 

The following parameters were considered: 


e 4 days (96 hours) injection period. 

¢ Flow rate at the wellbore was defined as 500 m®/day (5.79 x 10° m°/s). 

e Initial pressure was defined as 300 kef/ em”. 

e The wellbore radius was considered to be r,, = 0.108 m in all cases. 

e The vertical permeability was considered to be equal to the horizontal permeability for all cases. 
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e Porosity was defined as ¢ = 0.2 for all cases 


e The thickness of the shale between the layers (Ah) is considered to be zero in all cases for both the 
analytical and the numerical models. 


The other parameters can be found in table 1 for all cases considered: 

























































































Case | Properties Layer 1 Layer 2 Layer 3 
r(m) Region 1 Region 2 Region 3 | Region 1 Region 2 Region 3 Region 1 Region 2 | Region 3 
16 ee) 16 ee) = = - 
A k(mD) 1000 500 1000 500 - - 7 
h(m) 15 15 - 
bMo(cP) 5.1 - 
: Artificial : : Artificial : : , : 
ie Region 1 Beaune Region 3 | Region 1 pecene Region 3 Region 1 Region 2 | Region 3 
10 208 oo 10 208 oo - - = 
B k(mD) 500 550 550 500 500 550 - 7 
h(m) 10 10 ; 
Mo(cP) | 3.5 - 
rm) Region 1 Region 2 Region 3 | Region 1 Region 2 Region 3 Region 1 Region 2 | Region 3 
16 104 oo 10 208 oo 21 253 oo 
C k(mD) 500 450 400 500 450 400 500 450 400 
h(m) 15 15 15 
fio(cP) | 3.5 








Table 1 — Analyzed Cases 


In the graphs of all cases, besides the pressure change curves for the analytical and numerical solution, 
the curves of pressure derivatives with respect to the logarithm of time is also present. It is important to 
analyze the behavior of the pressure derivative as well, in order to properly interpret the results of the test. 

Consider the curve composed by blue circles to be the analytical solution for pressure variation and 
the curve composed by red circles to be the analytical solution for pressure derivative. For the numerical 
solutions, yellow stars curve represent the pressure variations and purple stars the derivative curve. 

First, consider case A with two layers and two regions and the radii of the regions are equal. Pressure 
and pressure derivative curves for the analytical model solution and numerical solution are presented in 





























figure 5: 
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Fig. 5 — Analytical and numerical pressure variation and derivative solutions for case A 
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It is possible to see a close agreement between the analytical and numerical simulated curves. 


Also, this graph represents clearly the two distinct regions of permeability impacting the pressure and 
derivative curves. 


Region 1 has a permeability value of 1000 mD and the second region has a permeability value of half of 
that, 500 mD. This directly affects the pressure change and its derivative curves. Indeed, for the initial 
times the pressure curve only notices the presence of the first region of permeability. After that, it reflects 
the second region, doubling the value of the derivative curve. 


The case B considers two layers and two regions like in the previous case, however, for this case the radii 
of the regions are distinct. This problem can be treated as a three regions of equal radii case. 


Case B graphs are given in figure 6: 
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dApwf (analytical) 
Apwf (simulated) 
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Fig. 6 — Analytical and numerical pressure variation and derivative solutions for case B 


The region of permeability 500 mD for layer 1 has 10m, and in layer 2, it has 228m and the region 
of permeability 550 mD for layer 1 has 228m, and in layer 2, it has 10m. Hence, the created artificial 
region 2 has a permeability value of 550 mD for layer 1 and 500 mD for layer 2 and for that region the 
presence of formation crossflow will have greater impact than in the other regions which have same values 
of permeability for both layers. 


Notice that the numerical simulated curves reflects a greater presence of the formation crossflow, since 
its curves are slightly below than in the analytical model’s. However, the behavior of the curves are still 
very similar. 


Finally, case C' is considered, three layers and three regions with distinct radii are considered. Since 
the radii for all regions are different, the number of different semipermeability coefficients is five for this 
example. 
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Case C’ graphs are given in figure 7: 
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Fig. 7 — Analytical and numerical pressure variation and derivative solutions for case C’ 


The numerical and analytical curves are relatively close. The fact that the value of permeability at 
the final regions for all layers is lower than the value at the initial regions causes the derivative curves to 
increase with time. 


In all three cases, it is possible to observe the proximity between the analytical solution and the numerical 
data, hence, the presented analytical model is able to describe the pressure behavior in multilayered radial 
composite reservoirs during the production periods. 


A. Equivalent Permeability. The pressure measured at the well is directly linked to the equivalent perme- 
ability (Keq). In that sense, another way to evaluating the efficiency of the model presented in this work is 
by estimating the equivalent permeability. 

For a multilayered reservoir, like presented in Cobb et al. (1972), it is given by: 


; fe) 

Keq = a [40] 
Do hy 
j=1 


n 


Since the model presented here reduces the cases of different regions radii into one with regions of equal 
radii, then the derivative profile enables the computation of equivalent permeabilities that combine the 
respective regions of the layers. That is, to obtain the first equivalent permeability, region 1 of layer 1 is 
combined with region 1 of layer 2 and so on. 

Using the source line logarithmic approximation, the following equation is used to estimate the equivalent 
permeability (Cobb et al., 1972): 





ae dLo 
Kea = 41 
1 Qhrm a] 


where m is the constant derivative level and hr is the total thickness. 
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In table 2, the estimated and real equivalent permeabilities are presented for the cases considered in the 
previous section along with the percentage errors: 


























Case a Keg real(mD) | Keg estimated (mD) Error (%) 
A Region 1 1000 1004.28 0.43 
B Region 1 500 504.37 0.87 
Region 2 525 524.05 0.18 
C Region 1 500 502.33 0.47 
Region 2 450 450.93 0.21 





Table 2 — Real and estimated equivalent permeability values and percentage error for cases A to C 


Most cases presented errors less than 1%. This indicates that the proposed formulation may be useful in 
obtaining reservoir parameters. 


B. Impact of Formation Crossflow. In order to observe the impact of formation crossflow in the pressure 
response, a case considering formation crossflow will be compared to one with the same reservoir properties 
however, without formation crossflow. For these cases, the vertical and horizontal permeabilities will be 


considered to be different. The values of the parameters can be found in Table 3: 
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Table 3 — Impact of Crossflow: Analyzed Cases 


To better perceive the impact of formation crossflow, very distinct values of permeabilities from one layer 
to another are considered. Figure 8 shows the graph results for the pressure change and pressure derivative 


for cases D, and Dg. 
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Fig. 8 — Pressure variation and derivative solutions for cases D, and Dp 
It is possible to see that after some time the curves become more apart from each other, that behavior 


can be explained by the impact of the formation crossflow. In figure 9, the pressure derivative curves can 
be observed closer: 
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Fig. 9 — Pressure derivative solutions for cases D, and D, 


The difference between the curves for cases D; and D2 are more noticeable in figure 9. 
Finally, table 4 presents the estimated and real equivalent permeabilities considering the first region of 
permeability, for cases D; and Dg, along with the percentage errors: 


Ke, estimated (nD) 


a 


Table 4 — Real and estimated equivalent permeability values and percentage error for cases D, and D2 





Case D,, which considers the presence of formation crossflow, has an estimated keg value closer to the 
real keg value. 
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5. Conclusions 


Based on the analytical solution proposed by Ehlig-Economides and Joseph (1987), this work proposed 
a formulation that combined the presence of formation crossflow with having different coupling regions of 
permeability along each layer for multilayered reservoirs under single phase flow. 

The model here suggested can be applied on a variety of cases. Cases with equal and different radii of 
regions of permeability, cases including a damaged region near the wellbore and it even can be extended to 
model cases including multilayered reservoir systems under two phase fluid flow. In that sense, the model 
here presented is quite robust and dynamic. 

For a set of cases, a comparison was made between the analytical solution and numerical simulation, and 
it was possible to see the agreement between them. The pressure response, along with other features, such 
as equivalent permeability, show a good agreement in the comparison for all cases. 


Nomenclature 


Ai, Bt - Coefficients for jth layer and ith region defined in the pressure solution 
ct - Total system compressibility 

h; - Thickness of layer j (m) 

hr - Total reservoir thickness (m) 

I, K; - Modified Bessel functions of first and second kind and order i 
keq - Reservoir equivalent permeability (mD) 

ks, - Horizontal permeability in layer j and region i (mD) 

kis - Vertical permeability in layer j and region i (mD) 

n - Number of layers in reservoir system 

m, - Number of regions in layer j 

p; - Initial Pressure (kgf/cm?) 

P; - Reservoir pressure in layer j and region i (kgf/cm?) 

Ap - Pressure change (kgf/cm?) 

p - Pressure change in the Laplace domain 

Pwf - Well-bottom hole pressure (kgf/cm?) 

q - Surface production rate (m?/day) 

qji - Flow rate for layer j and region i (m?/day) 

r - Radius (m) 

ri - Radius of region i in layer j (m) 

re - Reservoir outer radius (m) 

Tw - Wellbore radius (m) 

t - Time (h) 

u - Laplace variable 

X} - Semipermeability coefficient between layers 7 and j + 1 in region 7 
xi - Vertical permeability-thickness ratio 

K - Permeability-thickness product in layer j and region i 

oj; - Porosity in layer j 

ji; - Oil viscosity in layer j (cP) 

w; - Porosity-thickness product in layer j 
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